An explorative metabolomic analysis of the endothelium in pulmonary hypertension

Pulmonary hypertension (PH) is classified into five clinical diagnostic groups, including group 1 [idiopathic pulmonary arterial hypertension (IPAH) and connective tissue disease-associated PAH (CTD-aPAH)] and group 4 (chronic thromboembolic pulmonary hypertension (CTEPH)). PH is a progressive, life-threatening, incurable disease. The pathological mechanisms underlying PH remain elusive; recent evidence has revealed that abnormal metabolic activities in the endothelium may play a crucial role. This research introduces a novel approach for studying PH endothelial function, building on the genome-scale metabolic reconstruction of the endothelial cell (EC) to investigate intracellular metabolism. We demonstrate that the intracellular metabolic activities of ECs in PH patients cluster into four phenotypes independent of the PH diagnosis. Notably, the disease severity differs significantly between the metabolic phenotypes, suggesting their clinical relevance. The significant metabolic differences between the PH phenotypes indicate that they may require different therapeutic interventions. In addition, diagnostic capabilities enabling their identification is warranted to investigate whether this opens a novel avenue of precision medicine.

Pulmonary hypertension (PH) is a progressive disease with no available curative therapy, and despite recent advances in pharmacological therapy, the 5-year survival remains as low as 65-70 1,2 . Plasma metabolomic studies in idiopathic PAH (IPAH) have demonstrated that alterations in glucose homeostasis, lipid metabolism, and altered bioenergetics are associated with worse outcomes 3,4 . Metabolomic analysis of bone morphogenetic protein receptor type 2 mutations in the pulmonary endothelium revealed widespread metabolic reprogramming 5 , and recently, a phase 2 trial demonstrated that sotatercept, which targets dysfunctional bone morphogenetic protein pathway signalling, was efficacious in treating PAH 6 . The plasma metabolomic profile in chronic thromboembolic pulmonary hypertension (CTEPH) suggests aberrant lipid metabolism characterised by increased lipolysis, fatty acid oxidation, and ketogenesis 7 .
The endothelium is under normal conditions in a quiescent state. When activated, the endothelium secretes different growth factors and cytokines that affect endothelial cell (EC) proliferation, apoptosis, and coagulation; attract inflammatory cells; and/or affect vasoactivity to restore homeostasis 8 . Chronic activation of the endothelium leads to EC dysfunction, resulting in pathological changes, such as those observed in PH 9,10 . Many factors have been suggested to trigger EC dysfunction in PH, such as shear stress, hypoxia, inflammation, cilia length, and genetics [11][12][13] . The overactivated endothelium in PH secretes vasoconstrictive 14,15 and proliferative factors [15][16][17] and reduces the secretion of vasodilators, which indicates that EC dysfunction might play a central role in the pathogenesis of PH.
This study introduces a novel approach for investigating intracellular metabolism using a genome-scale metabolic reconstruction of the microvascular EC 18 . This EC model is constructed from the metabolic models of humans (RECON1), metabolomics data from human umbilical vein endothelial cells, a manual literature review of endothelial metabolism, and finally, transcriptomic data from three subtypes of endothelial cells (human umbilical vein, human mammary vascular, and human pulmonary artery) making it relevant for investigating the PH population 19,20 . Intracellular metabolic activities analysis. A complete overview of the results from the examined cellular activities are appended in Supplement Material A. Nineteen cellular activities had constant values of either 0 flux or 1000 flux indicating inactive or fully activated pathway across the samples, and were removed from further analysis.
When investigating the PH patients' profiles on the basis of cellular metabolic activities, four distinct clusters (Phenotypes A-D) were revealed that were independent of clinical PH diagnosis ( Fig. 2A). PLS-DA analysis revealed that the top 20 metabolic cellular activities separating the four phenotypes were NADH generation (VIP score above 4.5) followed by the degradation of 19 amino acids ( Clinical characteristics of the metabolic cellular activities of phenotypes A-D. The clinical characteristics and catecholamine measurement of patients with each phenotype based on metabolic cellular activities activity are presented in Table 2. Phenotype C had a median fourfold higher level of NT-proBNP than phenotype A, and phenotype D had a 14-fold higher level of NT-proBNP than phenotype B; phenotypes C and D also had a higher proportion of patients with WHO functional classifications III and IV.

Comparison of cellular metabolic pathway activities in phenotypes A-D. Amino acid metabo-
lism. Compared with phenotypes A and B, phenotype D had significantly increased degradation of amino acids (alanine, glycine, serine, histidine, asparagine, phenylalanine, valine, proline, tyrosine, aspartate, threonine, isoleucine, tryptophan, glutamate, leucine, glutamine, cysteine, ornithine, and lysine), indicative of increased catabolism; phenotype C had significantly increased degradation of amino acids expect for ornithine and lysine compared with phenotype B (Supplementary Material C).
Energy metabolism. Phenotype B had the highest ATP generation from glucose, and phenotype C had lower ATP generation from glucose than phenotype A (Fig. 3A). Phenotype C had the highest and phenotype B had the lowest proportion of ATP generation from the tricarboxylic acid (TCA) cycle (Fig. 3B).
TCA metabolism. Phenotype B had the lowest NADH generation by the TCA cycle, whereas phenotypes C and D had the highest (Fig. 3C).
Lipid metabolism. Phenotype C had lower activity in malonyl-CoA synthesis than phenotypes A and B, and phenotype A had lower activity than phenotypes B and C (Fig. 3D).

Discussion
The main finding of the present investigation was that the intracellular metabolic activity of ECs of PH patients clustered into four phenotypes (A-D) that were independent of the PH diagnosis. Notably, the disease severity (reflected by NT-proBNP levels) differed significantly between the metabolic phenotypes, suggesting their clinical relevance 21 ; phenotype D had the highest levels, and phenotype B had the lowest levels.
Our observations are in alignment with those of Kariotis and colleagues, who identified three whole-blood transcriptomic patient subgroups of IPAH patients that accounted for more than 90% of the cohorts 22 . These subgroups were associated with poor, moderate, and good prognoses.. Our study also extends to CTEPH patients, suggesting that the metabolic phenotypes are not restricted to PAH. Because the present study is, to the best of our knowledge, the first to explore the inferred intracellular endothelial metabolism in PH patients in the context of a genome-scale metabolic model (GEM), making direct comparisons with the literature is difficult. EC research in PH largely builds on cell culture experiments, which use ECs from PH patients undergoing lung transplantation or from biorepositories 23 . Unfortunately, these studies have been conducted under static conditions without hemodynamic shear stress, and the manipulated cells have been cultured on plastic, which is five orders of magnitude stiffer than in vivo conditions 23 ; thus, a direct comparison with results from metabolic systems analysis is difficult. Similarly, PH animal models poorly mimic complex human disease 24 and require extensive manipulation of the cells. This may explain the lack of previous reports on the identification of clinically relevant EC phenotypes with distinctly different intracellular metabolic activities.
In our exploratory modelling analyses, phenotype B had the lowest overall ATP generation but displayed the highest ATP generation from glycolysis. Prioritising glycolysis over oxidative phosphorylation may enable ECs to increase lactate production and thus promote angiogenesis [25][26][27][28] , which has previously been identified as a vital feature in studies of general metabolism in PAH vascular cells 29,30 . Utilising the glycolytic path may also minimise reactive oxygen species (ROS) production produced by oxidative phosphorylation in the TCA cycle and generate ATP more rapidly than the TCA cycle. Similar metabolic reprogramming has also been reported in cancer, providing ECs with a survival advantage 31 . The finding that phenotype B had the lowest NT-proBNP level and thus the mildest disease severity in the PH patients included in this study suggests that these ECs are subject to the least cellular stress.
Phenotype D had the highest disease severity as measured by NT-proBNP, which may be a consequence of high cellular stress, and thereby, elevated intracellular reactive oxygen species force the cell to increase amino acid degradation 32 . This is consistent with our findings in this study that amino acid degradation, particularly increased in phenotype D, was 19 of the top 20 pathways to distinct the phenotypes. Also, phenotype D had the highest levels of circulating catecholamines. In alignment with this, neurohumoral signalling, including sympathetic activation, was recently reported to contribute to the PAH pathophenotype 33 . Furthermore, Nagaya and www.nature.com/scientificreports/ colleagues investigated 60 PH patients and reported that plasma norepinephrine was significantly increased in functional class IV compared with functional classes II and III 34 . Conversely, in phenotype A, which exhibited severalfold lower NT-proBNP levels than phenotype D, significantly lower amino acid degradation and higher β-oxidation were observed. In addition, phenotype C, which had the second-highest disease severity as measured by NT-proBNP levels, displayed inverse ATP generation compared with phenotype B, with the lowest ATP generation from glycolysis and the highest oxidative phosphorylation, providing fuel to the TCA cycle from high amino acid degradation and beta-oxidation of fatty acids induced by low synthesis of malonyl-CoA.. Whether the differences in intracellular metabolic activities of ECs between the different PH phenotypes reflect differences in disease severities or genetic variation remain to be investigated.
The present study has several limitations inherent to its observational design; thus, no causal effect can be inferred. In addition, a limited number of PH patients from a single centre were studied, and this may further limit the generalisability of the results. Importantly, the PH patients received various pharmacological therapies with agents potentially affecting the endothelium, and the observed metabolic perturbations may be influenced by this. However, sample sizes will have to be considered, and validation of treatment impact on phenotype will await larger multicenter studies. Although the iEC3006 is the most comprehensive metabolic reconstruction of the EC to date, expansion of the metabolic coverage focusing on membrane lipids and glycan metabolism is required. Furthermore, the exploratory modeling analyses were based on a limited number of extracellular metabolites in the plasma, which may affect the intracellular metabolic cellular activity analysis despite covering all parts of the central metabolism 35 . Consequently, the inferred intracellular pathways analyses need to be validated in vitro under physiological settings. Also, metabolites from other cell types in plasma cannot, at the moment, be incorporated in the iEC3006 since it is a genome-scale single-cell metabolic model. However, our approach that plasma metabolism primarily reflects the EC alteration is based on the fact that 1 trillion EC is in constant contact with the circulating blood. Lastly, since SNP data were unavailable, we could not evaluate if the phenotype may reflect a continuum of PH progression severity or genetic differences.

Conclusions
The present exploratory study suggests that the endothelium displays four distinctly different intracellular metabolic phenotypes independent of the patient's PH diagnosis. The phenotypes are characterised by differences in ATP generation associated with varying disease severity as evaluated by NT-proBNP. The substantial metabolic differences between the PH phenotypes suggest that they may require different therapeutic interventions; furthermore, diagnostic capabilities enabling their identification is warranted to investigate whether this opens a novel avenue of precision medicine 36 . Patient selection. Patients were randomly selected retrospectively for enrollment in this study on the basis of their PH diagnosis from a biorepository of over 450 patients. All patients had the PH diagnosis confirmed with right heart catheterizations, including measurements of pulmonary pressures and cardiac output/index using the thermodilution method. In total, 12 patients with IPAH, 12 patients with CTD-aPAH, and 11 patients with CTEPH were selected, along with 20 healthy volunteers that have previously been published 18 . Clinical data from the patients were retrieved from a research database and electronic health records.

Methods
Analysis of clinical characteristics. Statistical analysis was performed using Rstudio (version 3.6.3).
Patients' descriptive data are presented as medians with interquartile ranges (IQRs) or as percentages (%). The nonparametric Kruskal-Willis test was used to evaluate differences in the main cellular metabolic activities in the four metabolic phenotypes.
NT-proBNP was adjusted to the standard maximum reference value as a percentage (i.e. the maximum reference value was equal to 100%).

Sample Preparation.
Blood samples were collected in EDTA tubes. Tubes were immediately centrifugated at 1,800 × g for 10 min at 5 °C two times to separate plasma. Plasma was aliquoted and frozen at − 80 °C.
Mass spectrometry analysis. Ultra-high performance liquid chromatography-mass spectrometry (UHPLC-MS) analysis was performed as previously published 37 using the Vanquish UHPLC system (Thermo Fisher Scientific, San Jose, CA, USA) with a Q Exactive mass spectrometer (HF Hybrid Quadrupole-Orbitrap, Thermo Fisher Scientific, San Jose, CA, USA). An electrospray ionisation interface was used as an ionisation source. The analysis was performed in negative and positive ionisation modes. A QC sample was analysed in www.nature.com/scientificreports/ MS/MS mode for the identification of compounds. UPLC was performed using the protocol described by Catalin et al., slightly modified so that the derivatisation reaction was stopped by the addition of chloroform 38 . Peak areas were extracted using Compound Discoverer 2.0 (ThermoFischer Scientific, Waltham, MA, USA). Gas chromatography-mass spectrometry (GC-MS) was used to detect amino and non-amino organic acids (7890B, Agilent, Santa Clara, CA, USA) coupled with a quadrupole mass spectrometry detector (5977B, Agient). The system was controlled by ChemStation (Agilent). Raw data were converted to netCDF format using Chemstation (Agilent) before the data was imported and processed in Matlab R2018b (Mathworks, Natick, MA, USA) using the PARADISe software 39 . Analysis of mass spectrometry data. On the basis of our prior work in critically ill patients and the ability to integrate the metabolites into the iEC3006, we selected 51 metabolites to be quantified (Supplementary Table D) 18 .
Less than 1% of the total values were missing, and these were imputed using the Missforest package 40 in RStudio (version 3.6.3). The quantified metabolic data were normalized by log10 transformation and Pareto scaling to create Principal Component Analysis (PCA), and the Volcano plot analysis. P-values were FDR adjusted.
Analysis of data with the iEC3006 genome-scale metabolic model (GEM). To date, our GEM EC, entitled iEC3006, is the most expansive genome-scale metabolic model of the EC, including 2,035 genes and 3,006 reactions involving a total of 2,114 metabolites 18 . Quantified metabolic data were analysed with iEC3006 using the Cobratoolbox in Matlab R2017b 41 . First, constraints on the uptake and secretion of metabolites were determined by the upper and lower quartiles of their respective transport reaction flux distributions as determined by random sampling flux analysis of the cultured human umbilical vein endothelial cell (HUVEC) constrained version of iEC3006.
We previously determined the normal metabolic variances by including 20 healthy volunteers variances in order to incorporate metabolic patient data into EC-GEM 18 . Mean fold changes (patient phenotype/healthy volunteers) from the plasma metabolomics data set were then applied to the upper and lower quantiles of each transport reaction to define the uptake and secretion in the patient phenotype models; if necessary, reactions were relaxed to obtain a feasible model (Supplementary Table E and Supplementary Table F). In total, 190 manually curated central metabolic cellular activities were investigated (Supplementary Table G) 35 using the metabolomics package in Rstudio (version 3.6.3). Constant values of metabolic cellular activities measurement were removed from further analysis. The intracellular metabolic similarity was identified by applying a heatmap with a dendrogram. The Euclidian distance and the complete cluster algorithm were used.

Data availability
Code to reproduce modeling analysis is available in the GitHub repository https:// github. com/ HHEN0 042/ PH